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SUMMARY 

Calculations of unsteady flows using a simplified marker and cell 
(SMAC), a pressure- implicit splitting of operators (PISO), and an 
iterative time -advancing scheme (ITA) are presented. A partial 
differential equation for incremental pressure is used In each 
time -advancing scheme. Example flows considered are a polar cavity flow 
starting from rest and self-sustained oscillatory flows over a circular 
and a square cylinder. For a large time- step size, the SMAC and ITA are 
more strongly convergent and yield more accurate results than PISO. The 
SMAC is the most efficient computationally. For a small time-step size, 
the three time -advancing schemes yield equally accurate Strouhal numbers. 
The capability of each time -advancing scheme to accurately resolve 
unsteady flows Is attributed to the use of a new pressure correction 
algorithm that can strongly enforce the conservation of mass. The 
numerical results show that the low frequency of the vortex shedding Is 
caused by the growth time of each vortex shed into the wake region. 


*Resident Research Associate at NASA Lewis Research Center. 
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coefficient for velocity of discrete momentum equation 
side of square cylinder 

coefficient for temporal discretization (i-1,2,3) 
frequency, f-l/T 
reference length 
pressure 

radius of circular cylinder 

normal distance measured from bottom wall of polar cavity 

period of vortex shedding 

time 

reference velocity 
velocity component, u^“{u,v) 
cartesian coordinates 
increment 

angle (unit in degree) 
molecular viscosity 
density 
summation 

non-dimensional time, r - tU 0 /L 
time -level 

neighboring grid points 
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INTRODUCTION 


For incompressible flows, the conservation of mass acts as a 
constraint condition that the velocity field needs to satisfy; while in 
compressible flows, the conservation of mass is given as a partial 
differential equation for the temporal variation of density. Due to this 
distinct difference, recent developments of numerical methods to solve 
compressible flows are mostly concentrated on unsteady Navier- Stokes 
equations while those for incompressible flows are mostly concentrated on 
steady Navier-Stokes equations. Thus the numerical methods to solve 
unsteady incompressible flows have been reported only sporadically even 
though one of the earliest numerical methods to solve an unsteady, 
incompressible flow appeared as early as 1965 [1] . 

A careful examination of various unsteady incompressible flow 
solution techniques can provide a valuable guideline for a further 
extension of these numerical methods to solve more complex unsteady flow 
problems. The accuracy, convergence nature, and computational effort of 
SMAC [2], PISO [3,4] and ITA are examined by solving a polar cavity flow 
[5] starting from rest and self - sustained oscillatory flows over a 
circular cylinder [2] and a square cylinder [6]. 

The ITA is a direct extension of steady flow solution techniques to 
solve unsteady flows. In the method, a steady flow solution technique is 
used iteratively to obtain a converged solution for each time -level. 
Therefore, the ITA can be computationally heavy. 

The marker and cell (MAC) scheme was originally developed to solve 
an unsteady free -surface flow. The method was simplified to solve an 
unsteady laminar flow over a circular cylinder by Braza et al. [2] and a 
few other flow cases cited in the reference. The name "marker and cell" 
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has been retained in Ref. [2] and in the present study even though the 
numerical method itself has nothing to do with the "marker and cell" for 
the unsteady flows considered in these studies. Two slightly different 
schemes, varying in the complexity of the pressure equation, have been 
presented in Ref. [2], The more complex case requires almost the same 
computational effort as that of PISO while the accuracy remains the same 
as that of the simpler case, and hence only the simpler case is 
considered herein. The simpler case is recast in a discrete form so that 
the context may be consistent with the other time -advancing schemes. 

The PISO [3] was proposed to avoid the heavy computational effort of 
ITA. Calculations of a laminar flow through a suddenly expanding pipe and 
a compressible flow entering a suddenly expanding closed pipe can be 
found in [4]. In the present study, the PISO [3,4] is modified to be 
memory efficient, which results in a slightly increased computational 
effort than the original PISO, For clarity, the modified PISO is 
described in detail in this paper. 

In the pressure correction algorithm used in this paper, the 
incremental pressure is obtained by solving a partial differential 
equation for incremental pressure, which yields a strongly diagonally 
dominant system of equations for the incremental pressure. It has been 
shown In a number of numerical calculations of various two- and 
three-dimensional flows that the present method yields strongly and 
highly convergent results even when highly graded and skewed meshes are 
used to discretize the flow domain [7-10], For example, the pressure 
correction algorithm yields a grid independent solution for a 3-D curved 
duct flow with a very small number of grid points and it can reduce the 
mass imbalance at least a few orders of magnitude smaller than those 
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obtainable using various other pressure correction algorithms [9]. 


NUMERICAL METHODS 

The incompressible laminar flow equations are given as; 


d\ii 


8xa 


- 0. 


d(pu^) d d 

+ (pUjUj) - 


at 


ax^ 


ax 4 


au^ au^ 

M( + — ) 

dxx dx L J 


> - — 


3p 

axi 


(i) 


( 2 ) 


where pu^u^ and pu^uj (i^j) represent the nonlinearity in each component 
of the momentum equation and the nonlinear coupling between the u^- and 
uj -velocity, respectively. Repeated indices imply summation over the 
index unless otherwise stated. 

The unsteady flow solution techniques are implemented on a 
pressure -staggered mesh and the incremental pressure is obtained by 
solving a partial differential equation for incremental pressure. Only a 
few necessary details of the present pressure correction algorithm are 
summarized in "Iterative Time -Advancing Scheme" sub-section even though 
the same pressure correction algorithm is used for each of the unsteady 
flow solution techniques. Details on the pressure correction algorithm 
can be found in Refs. [7-9]. 

The transient term in the momentum equation is treated implicitly. 
The temporal variation of pu^ can be discretized as 


d(P u i) n n-1 n-2 

- pC^u^ - p&2 u i. + P^3 u i 

at 


(3) 


5 


where (C^,C 2 ,C 3 ) -(1/At , 1/At , 0) and ( 3/2At , 4/2At , l/2At) for the first and 
the second order difference approximations, respectively. Since a few 
different time -advancing techniques are discussed below, the consistency 
of a few notations related to each time -advancing scheme is confined to 
each time -advancing scheme and each of such notations is explained only 
when it appears for the first time in each sub -section. 

Iterative Time- Advancing Scheme (ITA) 

• * ^ 

The initial guesses for the new time- level velocity (u^ ) and 

pressure (p*) are set equal to those of the previous time -level. The 

discrete momentum equation based on the guessed flow variables can be 

written as 

* * 3p* 

n-1 n-2 

(pCi + A i )u i - £ (A k u k } + + pCgUjL - pC 3 u i , 

nb 

no sum on i, (4) 

where u^ is a predicted velocity, is the coefficient of the 

-velocity, is a source term originating from a skewed mesh, and 

the coefficients A^ and are evaluated using the initial guess. The 
predicted velocities are obtained by solving eq. (4). The predicted 
velocities are not necessarily divergence free, and hence the velocities 
are corrected to satisfy the conservation of mass. Let u^*** and p** are 
the corrected velocity and pressure that satisfy the conservation of 
mass. Then 

- Uj_ + u' 1( (5) 

P - P + P (6) 
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where p # is the incremental pressure. The discrete momentum equations 
which satisfy the conservation of mass can be written as 


a ** 

* *** * *** * n-1 n-2 

(pC;L + A i )u i - £ {AfcU k } + + pC 2 u i - pCjUi , 

nb 3x^ 

no sum on i, (7) 

Subtracting eq. (4) from eq. (7) yields; 


1 3p' 

u^' - — , no sum on i, (8) 

(pC x + A*) dXi 


where the pressure gradient is left in continuous form deliberately. It 
is discussed later that the velocity-pressure decoupling that occurs when 
various pressure correction algorithms are used for a pressure-staggered 
mesh is caused by using a discretized pressure gradient In eq. (8) in 
deriving the discrete pressure correction equation. In deriving eq. (8), 
the summation over the neighboring grid points are disregarded. Retaining 
the residual originating from the neighboring grid points cause the 
numerical results to depend on the under -relaxation parameter [9] . 
Applying the conservation of mass to eq. (5) yields 

** 


3u'j 3uj 

3xj 3xj 


(9) 


Substituting eq. (8) into eq. (9) yields a partial differential equation 
for the incremental pressure given as; 


d [ 1 

ax^ [(pCj+A*) 
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where the last term in eq . (10) represents the mass imbalance. The 

incremental pressure is obtained by solving eq. (10) and the velocities 
and pressure are corrected using eqs . (5) and (6), respectively. The 

above procedure needs to be iterated using u^ and p as new initial 
guesses until the differences between the old and the new corrected 
velocities and pressures, respectively, become negligibly small; and the 
converged solutions represent the flow field of the current time -level. 

A single iteration of the ITA can not account for the nonlinearity 
in the momentum equation. However, the flow equations are solved 
iteratively until the convergence criteria are met at each time- level, 
the nonlinearity in each component of the momentum equation and the 
nonlinear coupling of u-and v-velocity are fully accounted for in the 
method. 

In each of the time -advancing schemes discussed in this paper, the 
incremental pressure is obtained by solving a partial differential 
equation for the incremental pressure. As all the central-differenced 
finite volume equations for self-adjoint second order elliptic partial 
differential equations are strongly diagonally dominant, the present 
discrete pressure correction equation is strongly diagonally dominant 
even for a highly skewed mesh. On the other hand, consider deriving a 
discrete pressure correction equation from eq. (8) using a discretized 
pressure gradient in eq. (8). Substituting eq. (8) into eq. (9) and 
integrating it over a pressure control volume yields a system of 
equations for incremental pressure which is not diagonally dominant. In 
such a case, the mass imbalance for a particular pressure grid point 
produces large pressure corrections for the adjacent pressure grid 
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points, and the velocity-pressure decoupling occurs. In case a discrete 
pressure gradient is used In eq. (8) to derive a system of discrete 
pressure correction equations, only a fully- staggered mesh or the 
momentum interpolation schemes [11-12] can yield a diagonally dominant 
system of equations for the Incremental pressure. However, the 
fully- staggered mesh is not an optimal grid layout for complex 
geometries; and the momentum interpolation scheme can not account for the 
grid skewness and it also yields a numerical result that depends on the 
under- relaxation parameter unless a specialized interpolation scheme is 
adopted [9,12]. However, the present pressure correction scheme does not 
yield a numerical result that depends on the under -relaxation parameter, 
since eq. (10) clearly states that the incremental pressure is driven 
only by the mass imbalance. 


Simplified Marker and Cell Scheme (SMAC) 

The SMAC and PISO are non-iterative time -advancing schemes and 
hence, an initial guess Is not necessary for these schemes. The discrete 
momentum equations based on the flow variables of the previous time -level 
can be written as 


n-1 * n-1 * n-1 

(pCi + )u£ - £ {A^ u^} + 

nb 


p„n- 1 

n 

+ pC2U£ 

dx i 


n-2 

pC 3 u i , 


no sum on i, (11) 


where u^* is the predicted velocity, A^ n ‘^ is the coefficient of the 
u^-velocity , S^ n '^ is a source term originating from a skewed mesh, and 
the coefficients Aj 11 " 1 and Sj 0 " 1 are evaluated using the flow variables 
of the previous time -level. The predicted velocities are obtained by 
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solving eq. (11). The predicted velocity field may not satisfy the 
conservation of mass, thus the predicted velocities need to be corrected 
to satisfy the conservation of mass. Let u^ n be the velocity that satisfy 
the conservation of mass. Then 

n * 

u i ’ u i - * , (12) 


where <f> is an auxiliary potential field. Taking divergence of eq. (12) 
yields 


d 

3xj 


d<j> 


3xj 


3uj 


3x4 


(13) 


which can be solved in the same way as that for eq. (10). The divergence 
free velocity field for the current time-level is obtained using eq. 
(12). In the SMAC, the search for the flow field of the current 
time- level is concluded by obtaining a consistent pressure (p ) that 
satisfies the momentum equation. The discrete momentum equation for the 
current time -level can be written as 


n- 1 n n- 1 n n-1 n-1 n-2 

(p )u^ - X u^} + S£ + pC — P^3 u i 

nb 5xj_ 

no sum on i 


(14) 


Subtracting eq. (11) from eq. (14) yields; 


n * 

pC 1 (u i - u t ) - 


d 

dx i 


( P n - p"- 1 ), 


( 15 ) 


where 


Ai n - 


has been disregarded since 


n-1 


«pC^ for a very small 
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time-step size. Substituting eq. (12) into (15) yields 


or 


3 d 

(pC]» (p n - p 11 ' 1 ) 


P n - P 11 ' 1 + PC\<f> 


( 16 ) 


where the spatial variation of has been neglected. 

In eq. (14), the coefficients of the U£ n -velocity and the source 
term originating from a skewed mesh have not been updated. Thus SMAC can 
not account for the nonlinearity in each component of the momentum 
equation. It also can not account for the nonlinear coupling of u- and 
v-velocity unless the discrete u- and v-momentum equations are solved 
simultaneously as a single system of equations. 

The velocity field obtained using the SMAC satisfy the conservation 
of mass strongly at each time-level. However, due to the use of an 
over-simplified pressure equation, eq. (16), the velocity and pressure 
fields may not satisfy the conservation of momentum rigorously. In fact, 
the over-simplified pressure equation can overshoot the pressure and It 
may yield a non-physical solution when the second order temporal 
discretization is used. On the other hand, the use of the first order 
temporal discretization always yields as accurate numerical results as 
those obtainable using the ITA or PISO incorporating the second order 
temporal discretization. This unusual behavior of the SMAC indicates that 
the second order temporal discretization is not quite compatible with the 
over-simplified momentum equation. 

A more elaborate SMAC presented in Ref. [2] has also been tested in 
the course of the present study, however, any significant Improvement in 
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accuracy was not observed. The same observation has also been expressed 
in Ref. [2]. On the other hand, the computational effort of the 
alternative scheme becomes comparable to that of the PISO. 

Pressure- Implicit Separation of Operators (PISO) 

The discrete momentum equation based on the flow variables of the 
previous time -level can be written as 


n-1 * * 5 P 

(pC± + )ui - H i - 


n-1 


n-1 n-2 

+ pC2U^ — P c 3 u i > no sum on i, (17) 




where 


* n-1 * n-1 

«i - Z (A k u k > + Si 


and the notations used In eq. (17) are the same as those for eq. (11). 

The predicted velocity is obtained by solving eq. (17). The predicted 
velocity field is not necessarily divergence free, and hence the velocity 
field needs to be corrected to satisfy the conservation of mass. Let u^ 
and p* are the first corrected velocity and pressure that satisfy the 
conservation of mass, respectively. Then 


** * 

* Uj_ + Uj/ , 


(18) 


P * - P «-i + p- 


(19) 


where u'i and p' are the first incremental velocity and pressure, 
respectively. The discrete momentum equations which satisfy the 
conservation of mass can be written as 


n-1 ** * dp n -l n -2 

{pCi + )u i - H i + pC 2 u i — pC 3 u t , no sum on i, (20) 

3x£ 


12 



Subtracting eq. (17) from eq. (20) yields; 


1 8p' 

Uj/ , no sura on i, (21) 

(pC^ + A^~l) ^ x i 

where the summation over the neighboring grid points and the source 
term originating from a skewed mesh have been disregarded in deriving 
eq. (21). Applying the conservation of mass to eq. (18) yields 

* 

<9u'.j du-j 

- - — ( 22 ) 

dxj <3xj 


Substituting eq. (21) into eq. (22) yields a partial differential 
equation for the incremental pressure given as; 


d 

dx 


j 


5p' 


[(pCi + A™" 1 ) dxA 
J J 


dui 


3x4 


(23) 


where the last term in eq. (23) represents the mass imbalance. The first 
incremental pressure is obtained by solving eq. (23), the first corrected 
pressure is obtained from eq. (19), and the first corrected velocity is 
obtained using eqs. (18) and (21). The first corrected velocity may 
satisfy the conservation of mass, but it does not satisfy the momentum 
equation according to Ref. [3,4], Let the velocity and pressure that 

tJLpiJLflJL* JL JL «JLi JL«JL «|U JU 

satisfy the momentum equation be u^ and p , where u^ and p are 
the second corrected velocity and pressure, respectively. The discrete 
momentum equation based on the second corrected velocity and pressure can 
be written as 
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** *** ** 
(pCl + A£ )u^ - 


+ pC 2 U £ 

3x| 


n- 2 

pC 3 U^ , no sum on i, (24) 


where 


** 

% 


** 

- I (A k u k } + 
nb 


** 
Si , 


and are evaluated using the first corrected velocity, and the 

momentum equation has been written in an explicit form to correct the 
velocity and pressure without solving a system of equations. Subtracting 
eq. (20) from eq. (24) yields; 

** , , ** * dp 

( pCi + )u^ - - Hj_ , no sum on i, (25) 

Sx^ 

where A^ 11 '^ can be set equal to A^** without incurring large error since 
Aj** « pC^ and 

f "kirk k"k 

u i “ u i - U£ ( 26 ) 

P" “ P** - P* (27) 


are the second incremental velocity and pressure, respectively. Inserting 
eq. (25) into eq. (26) and taking divergence of eq. (26) yields; 


d f 1 dp"' 

dxj [(pC^ + A**) dxj j 



(28) 


The second incremental pressure is obtained by solving eq . (28), the 
second corrected pressure is obtained from eq. (27), and the second 
corrected velocity is obtained from eqs. (25) and (26). 

The present PISO algorithm is slightly different from that of Issa 
[3] in the load vector term of eq. (28). In the original PISO algorithm, 
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the coefficients of the discrete momentum equation need to be stored to 
evaluate the load vector in eq. (29). In the present PISO algorithm, only 
the resultant given as needs to be stored to solve the second 

pressure correction equation. The present method is more memory efficient 
than the original PISO, while it requires more computational efforts. 

This modification has been motivated by the experience, obtained from 
numerical calculations of three-dimensional flows [9-10], that the 
required memory to solve complex three-dimensional flows may easily 
exceeds the current computer capacity. The modified PISO can better 
account for the nonlinearity in each component of the momentum equation 
since the updated residuals are used in solving the second pressure 
correction equation. 

In the first corrector step of the PISO, the velocity and pressure 
are corrected to satisfy the conservation of mass. In the second 
corrector step, the velocity and pressure are corrected to satisfy the 
momentum equation. Therefore, the driving force of the second pressure 
corrector step is obtained from the momentum imbalance even though the 
pressure correction equation is still derived from the conservation of 
mass. Consequently, the velocity field at each time-level may not satisfy 
the conservation of mass very accurately and this forms the fundamental 
difference between the SMAC and PISO. As in the SMAC, the PISO algorithms 
can not accurately resolve the nonlinear coupling of the u-and v-momentum 
equations unless the discrete u- and v-momentum equations are solved 
simultaneously as a single system of equations. 

NUMERICAL RESULTS 

Numerical results for the polar cavity flow [5] starting from rest 
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and self -sustained unsteady flows over a circular cylinder [2] and a 
square cylinder [6] obtained using the ITA, SMAC and PISO are presented 
below. Calculations of the example flows using the ITA and PISO show that 
the second order implicit time- stepping scheme yields stable numerical 
results for the time -step size as large as ten times of that for the 
first order implicit time-stepping scheme. The numerical results 
presented in this section are obtained using the first order temporal 
discretization for the SMAC and the second order temporal discretization 
for the ITA and PISO. 

The computational effort for each unsteady flow solution technique 
depends only slightly on each particular flow problem to be solved, the 
number of grid points and the number of time-steps, but it depends very 
strongly on the convergence criteria used. Since the system of equations 
are solved iteratively using a TDMA [13], one set of convergence criteria 
needs to be prescribed for the TDMA sweeps. The error norm used for the 
TDMA sweep is given as; 


1 

Mi 


l 

I R k > i/ R o , i I < C Z 


(29) 


where i - {u,v, or p) denotes each velocity component and pressure, R Q 
and are the sums of the absolute residuals of the discrete equation 
for every grid point evaluated at the Initial and at the k-th sweep of 
TDMA, and is the convergence criterion for the i-th flow variable. 
The SMAC and PISO are noniteartive time -advancing schemes, thus only the 
above set of convergence criteria needs to be prescribed. For the ITA, 
another set of convergence criteria is needed for the iterative solution 
of the flow equations. The error norm used for ITA is given as 
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( 30 ) 


\e\] - MAX {ABSKA^ - p jO/A^ i] j I j-1 ,N) < c \ 

where A n c £ denotes the maximum i-th flow variable at the n-th 
time -level, A n ^. j> is the ,0-th flow variable at the k-th iteration, and N 
denotes the number of grid points. In the ITA, the flow variables at each 
time- level evolve iteratively, thus an accurate solution for each 
discrete system of equations at each iteration is not necessary. Hence a 
much smaller maximum number of TDMA- sweeps than those for the SMAC and 
PISO can be assigned to avoid excessive TDMA- sweeps that can be caused by 
eq. (29). The convergence criteria used are cjp~ - {1x10"^, 1x10’^, 
lxl0‘2} and - {1x10*^, 1x10“ 1x10'^}. The numbers of maximum 
TDMA-sweeps for the SMAC and PISO are 11, 11, 100 for the discrete u- , 
v- , and p- equation, respectively; and those for ITA are 5, 5, 11 for the 
u- , v- , and p-equation, respectively. The maximum number of iterations 
for the ITA at each time -level is 11. The computational efforts (CRAY/YMP 
cpu-time) per each grid point and per each time step for the above 
convergence criteria are shown in Table 1. The use of slightly larger 
convergence criteria, smaller number of TDMA-sweeps and smaller number of 
iterations (for ITA) can decrease the computational effort significantly 
without losing the accuracy noticeably. 

Polar Cavity Flow Starting From Rest 

The lid-driven polar cavity flow is schematically shown in Fig. 
l-(a). The Reynolds number based on the lid velocity and the depth of the 
cavity is 350. The measured steady state velocity profiles can be found 
in Ref. [5]. The flow domain is discretized by 81x81 grid points. The 
evolution of u^ at (ry,9) - (0.246, 0.0) is shown in Fig. l-(b). The 


17 




normalized time for the polar cavity flow is based on the lid velocity 
and the depth of the polar cavity. For a large time-step size (Ar=0.05), 
the SMAC yields a slightly unstable numerical result at r~ 1 and the PISO 
does not yield a convergent solution. The calculated steady state 
velocity profiles are compared with the measured data in Fig. i-(c). It 
is shown in the figure that the velocity profiles obtained using 
different time -advancing techniques collapse into a single line at each 
9- location and that the numerical results are in very good agreement with 
the measured data. 

The mass imbalance for the polar cavity flow, obtained using the 
large time-step for each time -advancing scheme, is shown in Fig. 2. The 
mass imbalance produced by the ITA is one order of magnitude smaller than 
that produced by the SMAC. This result indicate that the ITA can most 
strongly enforce the conservation of mass and that the use of an 
under-relaxation do not obscure the numerical results. The 
ever-increasing mass imbalance produced by the PISO is mostly caused by 
the second corrector step for which the driving force is the momentum 
imbalance. For a small time-step size, the three methods yield equally 
accurate numerical results. However, the trend of the mass imbalance 
produced by the three methods remains the same as that shown In Fig. 2. 

Vortex Shedding Behind a Circular Cylinder 

A laminar flow over a circular cylinder at Reynolds number 100 is 
considered below. The Reynolds number is based on the diameter of the 
cylinder and the free stream velocity. A survey of measured data and a 
set of numerical results obtained using the SMAC can be found in Ref. 

[2]. In Ref, [2], the flow equations transformed onto a cylindrical-polar 


18 



coordinates were solved using a fully staggered mesh. In the present 
calculations, the inlet boundary is located at 8r upstream of the 
circular cylinder, the exit boundary is located at llOr downstream of the 
circular cylinder, and the side boundaries are located at 8r away from 
the circular cylinder. The flow domain is discretized by 146x101 grid 
points in x- , and y-coordinate directions, respectively. A uniform flow 
is prescribed at the inlet boundary and the vanishing gradient boundary 
condition is used at the exit boundary. A vanishing normal gradient of 
the tangential velocity and a vanishing transverse velocity are 
prescribed at the side boundaries. The flow field is perturbed by 
twisting the cylinder at the beginning so that the calculated velocity 
field can initiate the self -sustained oscillatory motion [2]. 

The streaklines calculated using the PISO are shown in Fig. 3. The 
streaklines show that the modified PISO can cleanly resolve the vortex 
shedding behind the circular cylinder. 

The calculated Strouhal numbers are compared with the measured data 
as well as the numerical results by Braza et al. [2] and Eaton [14] in 
Table 2. The normalized time for flow over a circular cylinder is based 
on the free stream velocity and the radius of the cylinder. It can be 
seen in the table that the Strouhal number (for Reynolds number of 110) 
obtained using a penalty finite element method [14] is smaller than than 
the measured data. The Strouhal numbers obtained using the three 
time -advancing schemes, implemented on a pressure-staggered mesh and 
using the new pressure correction algorithm, are in very good agreement 
with the measured data. 

The velocity vectors and the instantaneous streamlines passing 
through the separation and/or reattachment locations for the flow over 
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the circular cylinder is shown in Fig. 4. Instantaneous streamlines for 
unsteady flows are considered to be not meaningful by many researchers as 
yet. However, a close examination of the streamlines shown in Fig. 4 
reveals that the mass flow rate across the streamlines are by far smaller 
than that along the streamlines. Thus the streamlines passing through the 
separation and the reattachment locations can be considered as an 
enclosure of the eddy attached to the cylinder. This interpretation of 
the instantaneous streamline is similar to that for steady flows. It can 
be seen in the figure that the attached vortex grows slowly until it 
becomes almost as large as the circular cylinder, and then the fully 
grown vortex is swept into the wake region by the nearby free stream. 
These figures clearly indicate that the low frequency (or the long 
period) of the vortex shedding is caused by the growth time of the 
attached vortex. It is also shown in the figure that each vortex shed 
into the wake region is accompanied by an alternating switching of the 
reattaching streamline with one of the two instantaneous streamlines 
passing through the separation locations. Undoubtedly, the instantaneous 
streamlines help to better understand the vortex shedding mechanism. It 
is also shown later in this section that the streamlines help to clearly 
identify the separation and reattachment locations. 

The streamline contours and a carpet plot of stream function at 
t - 0.14T are shown In Figs. 5- (a) and 5-(b), respectively. The saddle 
point shown in Fig. 5-(b) is caused by a closely spaced counter- ratating 
vortices, see Ref. [15] for more details. 

The evolution of stagnation location in time for flow over a 
circular cylinder is shown in Fig. 6, It can be seen in Fig. 6-(b) that 
the numerical result obtained using the PISO exhibits a strong dependence 
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on the time-step size. The amplitude and the phase difference obtained 
using the P1S0 with Ar—O.OS is in worst agreement with the other 
numerical results including the one obtained using the same PISO with 
Ar=0.01. As the time -step size is decreased, both numerical results 
obtained using the SMAC and PISO approach that obtained using ITA. These 
results indicate that the ITA yields the most accurate results due to its 
capability to resolve the nonlinearity in each component of the momentum 
equation and the nonlinear coupling of the u- and v-velocity and the 
capability to enforce the conservation of mass most strongly. The strong 
dependence of the PISO on the time-step size is caused by the second 
corrector step which cause the velocity field less divergence free and by 
the linearized momentum equation that yields a solution regardless of the 
conservation of mass is rigorously satisfied or not. 

The time -varying separation and reattachment locations are shown in 
Fig. 7-(a) and 7-(b), respectively. For the flow over the circular 
cylinder, the mesh and the time -step size are small enough to accurately 
resolve the vortices generated around the smooth cylinder . Thus , the 
numerical results obtained using the ITA and SMAC exhibit only a very 
small phase difference, while the PISO yields a numerical result that 
deviates most from the other numerical results. The reattachment location 
presented in Ref. [2] is, in fact, another separation point located at 
the symmetric lower part as shown in Fig. 4. 

The lift and drag forces are shown in Figs. 8-(a) and 8-(b), 
respectively. As in Fig. 7, the numerical results obtained using the ITA 
and SMAC are In very good agreement with each other, while the PISO not 
only over-predict the amplitude but also yields a significant amount of 
phase difference. The present numerical methods yield slightly larger 


21 


amplitudes for the lift and drag than those presented in Ref. [2]. 

Vortex Shedding Behind a Square Cylinder 

A laminar flow over a square cylinder at Reynolds number 190 is 
considered below. The Reynolds number is based on a side of the square 
cylinder (b) and the free stream velocity. The measured Strouhal number 
and a numerical result obtained using the ITA can be found in Ref. [6]. 

In Ref. [6], the flow domain was discretized using a fully staggered mesh 
and the convection terms were discretized using a formally third order 
accurate QUICKEST scheme. In the present calculations, the inlet boundary 
is located at 5.5b upstream of the square cylinder, the exit boundary is 
located at 35b downstream of the square cylinder, and the side boundaries 
are located at 5.5b away from the square cylinder. The flow domain is 
discretized by 131x101 grid points in x-, and y-coordinate directions, 
respectively. The boundary conditions and the Initial perturbations are 
the same as those used for the flow over a circular cylinder. For the 
flow over the square cylinder, the numerical methods yield the 
self-sustained oscillatory motion without the use of the initial 
perturbation. However, it takes a while for the fully oscillatory motion 
to be established. 

The streaklines calculated using the SMAC are shown in Fig. 9. The 
streaklines show that the SMAC can cleanly resolve the vortex shedding 
behind the square cylinder. 

The calculated Strouhal numbers are compared with the measured data 
as well as the numerical results obtained by Davis and Moore using the 
QUICKEST scheme [6] in Table 3. The normalized time for flow over the 
square cylinder is based on the free stream velocity and a side of the 
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square cylinder. It can be seen in the table that the present numerical 
results are in very good agreement with the measured data, while the 
numerical results obtained by Davis and Moore deviate further as the mesh 
is refined. The calculated Strouhal number, which deviates farther from 
the measured data as the mesh is refined, is caused by the QUICKEST 
scheme [6]. 

The velocity vectors and the streamlines for the flow over the 
square cylinder is shown in Fig. 10. As in the flow over the circular 
cylinder, the attached vortex grows slowly until it becomes almost as 
large as the square cylinder, and then the fully grown vortex is swept 
into the wake region by the nearby free stream'. Thus the vortex shedding 
mechanism and the cause for the low frequency vortex shedding are 
essentially the same as those for the flow over the circular cylinder. 

However, in the present case, the attached vortex is split into two parts 

by the sharp corner and hence, there almost always exists another pair of 
separation and reattachment locations than the flow over the circular 
cylinder . 

The evolution of stagnation location in time for flow over a square 
cylinder is shown in Fig. 11. It can be seen in Fig. 11- (b) that the 
numerical results obtained using different time -advancing schemes yield 
small phase differences, but all the numerical results exhibit neatly 
organized oscillatory motion. 

The time -varying lift coefficient for flow over the square cylinder 
is shown In Fig. 12. It can be seen in the figure that the amplitude 

obtained using the PISO with Ar*=0.05 is in worst agreement with the other 

numerical results including the one obtained using the same PISO with 
Ar=0.01. As the time -step size is decreased, the amplitudes obtained 
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using the SMAC and PISO approach that obtained using the ITA. Again, the 
undesirably strong dependence of the PISO on the time -step size is caused 
by the second corrector step which cause the velocity field to deviate 
further from a divergence free velocity field. 

The drag for flow over the square cylinder is shown in Fig. 13. For 
the flow over the square cylinder, the time-step size is not small enough 
to accurately resolve the initially perturbed flow field. Thus, the 
numerical results obtained using the different time -advancing schemes 
exhibit large phase differences. However, all the numerical methods, 
except the PISO, yield almost the same frequencies and the amplitudes. 

CONCLUSIONS AND DISCUSSION 

Calculations of a polar cavity flow starting from rest and 
self-sustained oscillatory flows over a circular and a square cylinders 
using the ITA (Iterative Time -Advancing Scheme), SMAC (Simplified Marker 
and Cell) and PISO (Pressure Implicit Splitting of Operators) are 
presented. 

The numerical results show that the SMAC is the most efficient 
computationally and yields accurate results. Calculations of the 
lid- driven polar cavity show that the SMAC is even compete tive with 
steady flow solvers to solve steady flows while the ITA, PISO and many 
other unsteady flow solvers are not [9,15], 

The ITA can account for the nonlinearity in each component of the 
momentum equation and the nonlinear coupling of the u- and v- velocity 
through the iterative solution of the flow equations. Thus the ITA yields 
the most accurate numerical results for a large time-step size. The SMAC 
and PISO are non-iterative time -advancing schemes and hence these methods 
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can only weakly account for the nonlinearity in the Navier-Stokes 
equations through the use of predictor and corrector steps. As the 
time-step size is decreased, the numerical results obtained using the 
SMAC and PISO approach those obtained using the ITA, which shows that the 
SMAC and PISO can accurately resolve the nonlinearity in the 
Navier-Stokes equations if a sufficiently small time-step is used. 

The SMAC and PISO are quite different in their nature. In the SMAC, 
the predicted velocity field is corrected to satisfy the conservation of 
mass and the conservation of momentum is achieved by obtaining a 
consistent pressure. In the PISO, a second corrector step is Introduced 
to correct the momentum imbalance. The numerical results show that the 
second corrected velocity field deviates farther from a divergence free 
velocity field and that the linearized momentum equation yields a 
solution regardless of the conservation of mass is rigorously satisfied 
or not. Thus the PISO exhibits an undesirable strong dependence on the 
time-step size and it is also weakly convergent. These numerical results 
indicate that accurately enforcing the conservation of mass is very 
important to enhance the convergence nature as well as to obtain accurate 
numerical results. The capability of each time -advancing scheme to 
accurately resolve the unsteady flows is largely attributed to the 
capability of the new pressure correction algorithm which can strongly 
enforce the conservation. 

Numerical results for flows over a circular and a square cylinders 
show that a small vortex attached to the cylinder grows very slowly in 
time until it becomes as large as the cylinder and then the fully grown 
vortex is shed into the wake region. Each vortex shedding is accompanied 
by an alternating switching of the reattaching streamline with one of the 
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two streamlines passing through the separation locations, and the 
alternating switching generates a small attached vortex. Each vortex 
shedding constitute a half cycle of the self -sustained oscillatory 
motion. The vortex leaving the cylinder is of the same size as that of 
the cylinder, and hence the low frequency is caused by the growth time of 
the attached vortex. 
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TABLE I 


Computational Efforts for Each Unsteady Flow Solution Method 
(CRAY/YMP cpu-time per each grid point and per each time step) 


ITA 

SMAC 

PISO 

1.7xl0- 4 

2.5xl0' 3 

4. 5x10' 3 


(unit in second) 


TABLE II 

Strouhal Numbers (S c -2fr/U 0 ) for Flow over a Circular Cylinder 



ITA 

SMAC 

PISO 

Ref. [2] 

Ref. [14] 

* 

rt 

A r 

0.05 

0.05 

0.01 

0.05 

0.01 

0.01 





St 

1 

0.158 

0.155 

0.157 

0.164 

0.160 

0.16 

0.147 

0.16 


TABLE III 

Strouhal Numbers (S^— fb/U 0 ) for Flow over a Square Cylinder 



ITA 

SMAC 

PISO 

Ref. [6] 

Exp' t 

At 

0.05 

0.05 

0.01 

0.05 

0.01 

0.05 a 

0 . 05 b 

— 

St 

0.133 

0.126 

0.13 

0.142 

0.135 

0.159 

0.165 

0.146 


(a: 41x40 mesh, b: 51x62 mesh) 
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ITA(At=0.05) — — - PISO(At=0.05) SMAC(At=0.05) 

Figure 2.— Mass imbalance. 
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Figure 3 — Streaklines for flow over a circular cylinder. 
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leous streamline contours passing through the separation 
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(b) Carpet plot. 

Figure 5. — Streamlines at t - 0.14T. 
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(a)PISO (Ar - 0.01). 


I 

100 


125 



- PISO(At= 0.05) • — - PISO(At=0.01) 
(b) Comparison of iTA, SMAC and PISO. 


Figure 6.— Stagnation location for flow over a circular cylinder. 
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(b) Drag coefficient. 

Figure 8.— Lift and drag for flow over a circular cy finder. 
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Figure 9.— Streaklines for flow over a square cylinder. 
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Figure 10 — Velocity vectors and instantaneous streamlines passing through the separation a 
reattachment locations. 


38 

















STAGNATION LOCATION 








(a)SMAC (Ar *0.01). 



ITA(At=0.05)— - - SMAC(At= 0.05) — — SMAC(Ax=0.0 1) 

PISO(Ax=0.05) PISO(At=0.01) 

(b) Comparison of 1TA, SMAC and PISO. 

Figure 12.— Lift for flow over a square cylinder. 
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DRAG COEFFICIENT 



T 


(a)SMAC (Ar -0.01). 



X 


1 ITA(At=0.05) — SMAC(At=0.01) PISO(At=0.01) 

(b) Comparison of ITA, SMAC and PISO. 

Figure 13 — Drag for flow over a square cylinder. 
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